# Replication script for the example and plots in
# Iacus, King, Porro (2018), "A Theory of Statistical Inference for Matching Methods in Causal Research"
#
# These scripts have been tested under
# R version 3.4.4 (2018-03-15)
# CEM package Version: 1.1.17
#
# Any issue with this script? Please contact: stefano.iacus@unimi.it

library(cem)
data(LL)
library(MatchIt)


# exact match on X
LL2 <- LL
LL2$age <- sprintf("%.15f", LL$age )
LL2$education <-  sprintf("%.15f", LL$education )
LL2$black <-  sprintf("%.15f", LL$black )
LL2$married <-  sprintf("%.15f", LL$married )
LL2$nodegree <-  sprintf("%.15f", LL$nodegree )
LL2$hispanic <-  sprintf("%.15f", LL$hispanic )
LL2$u74 <-  sprintf("%.15f", LL$u74 )
LL2$u75 <-  sprintf("%.15f", LL$u75 )
LL2$re74 <-  sprintf("%.15f", LL$re74 )
LL2$re75 <-  sprintf("%.15f", LL$re75 )


exactOnX <- cem("treated", data=LL2,drop="re78")
exactOnX

# exact match on the pscore
logit <- glm(treated~age + I(age^2) + education + I(education^2) + black +
+ hispanic + married + nodegree + re74  + I(re74^2) + re75 + I(re75^2) +
+ u74 + u75, data=LL, family="binomial")

fit <- sprintf("%.15f", logit$fitted)
tab <- table(fit)
which(tab>1) -> idx
sum(tab[idx])
tab2 <- table(fit,LL$treated)
idx2 <- which(tab2[,1]>0 & tab2[,2]>0)

tab2[idx2,]
colSums(tab2[idx2,])


# cal <- 0.001 ## used later to extract information




set.seed(123)

m.out <- matchit(formula=treated~age + I(age^2) + education + I(education^2) + black +
hispanic + married + nodegree + re74  + I(re74^2) + re75 + I(re75^2) +
u74 + u75, data=LL, method = "nearest", caliper=0.1)

myvars <- c("treated", "age","education","black","hispanic","married","nodegree","u74","u75","re74","re75")

mymatches <- list(data.frame(names(m.out$w)))
names(mymatches[[1]])[1] <- "id"
mymatches[[1]]$weight <- m.out$w
mymatches[[1]]$method <- "matchit psm"

m.out

mmat <- m.out$match.matrix
idx <- which(!is.na(mmat))
matched <- c(rownames(mmat)[idx], mmat[idx])
ps <- m.out$distance[matched]
LL[matched,] -> LL.psm
dim(LL.psm)

mmat.ps <- mmat[idx,1]

for(i in names(mmat.ps)){
    lab1 <- i
    lab2 <- as.character( mmat.ps[i] )
    cat(sprintf("\n%.5f, %.5f, %.6f\n", ps[lab1], ps[lab2], abs(ps[lab1]- ps[lab2])))
    print(rbind( LL[lab1, myvars],  LL[lab2, myvars] ))
}
   
mat <- cem("treated", data=LL,drop="re78", cut=list(age=c(17, 29.67, 42.33, 55),education=c(3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), black=3, married=3, hispanic=3, nodgree=3, u74=3, u75=3, re74=3, re75=3 ))
mat

mat.psm <- cem("treated", data=LL.psm,drop="re78",cut=mat$breaks)
mat.psm

# labour force stata
matLFS <- cem("treated", data=LL,drop=c("re78","re74","re75"), cut=list(age=c(16, 24.5, 54.001, 65),education=c(0, 6.001, 12.001, 16.01, 19.01), black=3, married=3, hispanic=3, nodgree=3, u74=3, u75=3 ))
matLFS

matLFS2 <- cem("treated", data=LL,drop=c("re78"), cut=list(age=c(16, 24.5, 54.001, 65),education=c(0, 6.001, 12.001, 16.01, 19.01), black=3, married=3, hispanic=3, nodgree=3, u74=3, u75=3 ))
matLFS2



#LL$black <- as.character(LL$black)  # in case you want to force extact matching
#LL$hispanic <- as.character(LL$hispanic)
#LL$married <- as.character(LL$married)
#LL$nodegree <- as.character(LL$nodegree)
#LL$u74 <- as.character(LL$u74)
#LL$u75 <- as.character(LL$u75)

set.seed(123)
m.out1 <- matchit(formula=treated~age + I(age^2) + education + I(education^2) + black +
hispanic + married + nodegree + re74  + I(re74^2) + re75 + I(re75^2) +
u74 + u75, data=LL, method = "nearest", caliper=0.001)

m.out2 <- matchit(formula=treated~age + I(age^2) + education + I(education^2) + black +
hispanic + married + nodegree + re74  + I(re74^2) + re75 + I(re75^2) +
u74 + u75, data=LL, method = "nearest", caliper=0.01)

m.out3 <- matchit(formula=treated~age + I(age^2) + education + I(education^2) + black +
hispanic + married + nodegree + re74  + I(re74^2) + re75 + I(re75^2) +
u74 + u75, data=LL, method = "nearest", caliper=0.1)

m.out4 <- matchit(formula=treated~age + I(age^2) + education + I(education^2) + black +
hispanic + married + nodegree + re74  + I(re74^2) + re75 + I(re75^2) +
u74 + u75, data=LL, method = "nearest")




mm1 <- data.frame(names(m.out1$w))
names(mm1)[1] <- "id"
mm1$weight <- m.out1$w
mm1$method <- "psm 0.001"

mm2 <- data.frame(names(m.out2$w))
names(mm2)[1] <- "id"
mm2$weight <- m.out2$w
mm2$method <- "psm 0.01"

mm3 <- data.frame(names(m.out3$w))
names(mm3)[1] <- "id"
mm3$weight <- m.out3$w
mm3$method <- "psm 0.1"


mm4 <- data.frame(names(m.out4$w))
names(mm4)[1] <- "id"
mm4$weight <- m.out4$w
mm4$method <- "psm no caliper"

mm5 <- data.frame(rownames(LL))
names(mm5)[1] <- "id"
mm5$weight <- matLFS$w
mm5$method <- "cem LFS"

mm6 <- data.frame(rownames(LL))
names(mm6)[1] <- "id"
mm6$weight <- matLFS2$w
mm6$method <- "cem LFS2"


mymatches <- list(mm1, mm2, mm3, mm4,mm5,mm6)

mymatches <- list(mm1, mm2, mm3) # without optimal CEM, which is too good


set.seed(123)
sp1 <- spacegraph("treated", data=LL,drop="re78",R=list(cem=0), M=2, other.matches=mymatches)
sp2 <- spacegraph("treated", data=LL,drop="re78",R=list(cem=500), M=2)
sp3 <- spacegraph("treated", data=LL,drop="re78",R=list(cem=0,psm=300), psmpoly=2, M=2)
sp4 <- spacegraph("treated", data=LL,drop="re78",R=list(cem=0,mdm=300),  M=2)
spA <- combine.spacegraphs(sp1,sp2)
spB <- combine.spacegraphs(spA,sp3)
spC <- combine.spacegraphs(spB,sp4)

# if(interactive())
#  plot(spC,N="all", balance.metric="mdiff",xlim=c(730,200),ylim=c(0,0.11),scale.var=T,main="")

save.image("all.rda")

space <- spC$space

pch.mdm <- 15
pch.psm <- 8
pch.cem <- 20
col.psm <- rgb(0,0,0,alpha=0.5) #"black" # rgb(0,1,0,alpha=0.5)
col.mdm <- rgb(0,0,0,alpha=0.5) #"black" # rgb(1,0,0,alpha=0.5)
col.cem <- rgb(0,0,0,alpha=0.5) #"black" # rgb(0,0,1,alpha=0.5)

idxsp.raw <- which(space$Method=="raw")
idxsp.psm <- which(space$Method=="psm")
idxsp.mdm <- which(space$Method=="mdm")
labs.psm <- c("psm 0.001", "psm 0.01", "psm 0.1")
idxsp.psm2 <- which(space$Method %in% labs.psm)

space$tot <- space$G0+space$G1
space$dotpch <- pch.cem
space$dotpch[idxsp.psm] <- pch.psm
space$dotpch[idxsp.mdm ] <- pch.mdm
space$dotpch[idxsp.raw] <- 21
space$dotpch[idxsp.psm2] <- 21

space$dotcol <- col.cem
space$dotcol[idxsp.psm] <- col.psm
space$dotcol[idxsp.mdm] <- col.mdm
space$dotcol[idxsp.raw] <- "black"


# FIGURE 1 in the paper, file fig1-bw.pdf (Color) in the paper

pdf("fig1-bw.pdf",width=11,height=5,pointsize=12)
# to produce png replace the revious line with
# png("fig1-bw.png",width=900,height=600,pointsize=12)
# the same for all other figures below

par(mar=c(4,4,1,1))
plot(1/(space$G0+space$G1),space$Mdiff,
pch=space$dotpch,
col=space$dotcol,
xlim=c(1/730,1/180),ylim=c(0,0.11),cex=0.5,
ylab = "Average Difference in Means",
xlab = "number matched, scaled as 1/matched",
axes=FALSE)
axis(2)
axis(1,at=1/space$tot,labels=space$tot)
box()

legend(1/(space$tot)[idxsp.raw],space$Mdiff[idxsp.raw], "raw", box.col="black",bg="white",adj=0,cex=.5,
x.intersp = 0.1,y.intersp = 0.1,seg.len=0)


for(i in 1:3)
 legend(1/(space$tot)[idxsp.psm2][i],space$Mdiff[idxsp.psm2][i], labs.psm[i],box.col="black",bg="white",adj=0,cex=.5,
x.intersp = 0.1,y.intersp = 0.1,seg.len=0)

# best stratification solutions [graphically inspected]
ii1 <- which(space$tot>630 & space$Mdiff < 0.01)
ii2 <- which(space$tot>680 & space$Mdiff < 0.015)
ii <- c(ii1,ii2)
points(1/space$tot[ii],space$Mdiff[ii])

text(1/720,0.005, labels="best\nrandom\nstratifications",cex=0.5)

dev.off()


# FIGURE 1 in the paper, file fig1-bw.pdf (Color) in the paper

pdf("fig1-col.pdf",width=11,height=5,pointsize=12)

col.psm <-   rgb(0,1,0,alpha=0.5)
col.mdm <-   rgb(1,0,0,alpha=0.5)
col.cem <-   rgb(0,0,1,alpha=0.5)

space$dotcol <- col.cem
space$dotcol[idxsp.psm] <- col.psm
space$dotcol[idxsp.mdm] <- col.mdm
space$dotcol[idxsp.raw] <- "black"

par(mar=c(4,4,1,1))
plot(1/(space$G0+space$G1),space$Mdiff,
pch=space$dotpch,
col=space$dotcol,
xlim=c(1/730,1/180),ylim=c(0,0.11),cex=0.5,
ylab = "Average Difference in Means",
xlab = "number matched, scaled as 1/matched",
axes=FALSE)
axis(2)
axis(1,at=1/space$tot,labels=space$tot)
box()

legend(1/(space$tot)[idxsp.raw],space$Mdiff[idxsp.raw], "raw", box.col="black",bg="white",adj=0,cex=.5,
x.intersp = 0.1,y.intersp = 0.1,seg.len=0)


for(i in 1:3)
legend(1/(space$tot)[idxsp.psm2][i],space$Mdiff[idxsp.psm2][i], labs.psm[i],box.col="black",bg="white",adj=0,cex=.5,
x.intersp = 0.1,y.intersp = 0.1,seg.len=0)

# best stratification solutions [graphically inspected]
ii1 <- which(space$tot>630 & space$Mdiff < 0.01)
ii2 <- which(space$tot>680 & space$Mdiff < 0.015)
ii <- c(ii1,ii2)
points(1/space$tot[ii],space$Mdiff[ii])

text(1/720,0.005, labels="best\nrandom\nstratifications",cex=0.5)

dev.off()



# FIGURE 1 in the paper, interactive
# if(interactive())
# plot(spC,N="all", balance.metric="mdiff",xlim=c(730,200),ylim=c(0,0.11),scale.var=TRUE,main="")


# ATT estimation for different matching solutions

obj <- list(data=LL, match=cem(treatment="treated",data=LL, drop="re78",eval.imbalance=FALSE))

myf <- formula(re78~treated+age+education+black+hispanic+married+nodegree+u74+u75+re74+re75)
# PSM
idx.psm <- which(spC$space$Method=="psm")

psm <- cem:::psm

all.psm <- NULL
all.psm.sd <- NULL
all.psm.n <- NULL

k <- 0
for(i in idx.psm){
    k <- k+1
    est <- eval(parse(text=spC$space$Call[i]))
    tmp <- lm(myf , data=est$match.dat, weights=est$match.dat$weights)

    temp.est <- coef(tmp)["treated"]
    temp.sd <- coef(summary(tmp))["treated","Std. Error"]
    temp.n <- spC$space$G0[i]+spC$space$G1[i]
    cat(sprintf("\nPSM - %.3d : %f (%f) n=%d", k,temp.est,temp.sd,temp.n))

    all.psm <- c(all.psm, temp.est)
    all.psm.sd <- c(all.psm.sd,temp.sd)
    all.psm.n <- c(all.psm.n, temp.n)


}
tab.psm <- data.frame(est=all.psm, sd=all.psm.sd, n=all.psm.n,method="psm")


save(tab.psm, file="psm.rda")
#load("psm.rda")

summary(all.psm)
if(interactive())
 plot(density(all.psm))

# MDM
idx.mdm <- which(spC$space$Method=="mdm")

mdm <- cem:::mdm

all.mdm <- NULL
all.mdm.sd <- NULL
all.mdm.n <- NULL

k <- 0
for(i in idx.mdm){
    k <- k + 1
    est <- eval(parse(text=spC$space$Call[i]))
    tmp <- lm(myf, data=est$match.dat, weights=est$match.dat$weights)

    temp.est <- coef(tmp)["treated"]
    temp.sd <- coef(summary(tmp))["treated","Std. Error"]
    temp.n <- spC$space$G0[i]+spC$space$G1[i]
    cat(sprintf("\nMDM - %.3d : %f (%f) n=%d", k,temp.est,temp.sd,temp.n))

    all.mdm <- c(all.mdm, temp.est)
    all.mdm.sd <- c(all.mdm.sd,temp.sd)
    all.mdm.n <- c(all.mdm.n, temp.n)

}
tab.mdm <- data.frame(est=all.mdm, sd=all.mdm.sd, n=all.mdm.n,method="mdm")

save(tab.mdm, file="mdm.rda")
#load("mdm.rda")

summary(all.mdm)
if(interactive())
 plot(density(all.mdm))


# CEM

idx.cem <- which(spC$space$Method=="cem")
all.cem <- NULL
all.cem.n <- NULL
all.cem.sd <- NULL
k <- 0

for(i in idx.cem){
    k <- k+1
    xx <- eval(parse(text=spC$space$Call[i]))
    est <- att(xx, re78~treated+age+education+black+hispanic+married+nodegree+u74+u75+re74+re75, data=LL)
    temp.est <- est$att.model["Estimate","treated"]
    temp.sd <- est$att.model["Std. Error","treated"]
    temp.n <- spC$space$G0[i]+spC$space$G1[i]
    cat(sprintf("\nCEM - %.3d : %f (%f) n=%d", k,temp.est, temp.sd, temp.n))
    all.cem <- c(all.cem, temp.est)
    all.cem.sd <- c(all.cem.sd, temp.sd)
    all.cem.n <- c(all.cem.n, temp.n)
}

tab.cem <- data.frame(est=all.cem, sd=all.cem.sd, n=all.cem.n,method="cem")

save(tab.cem, file="cem.rda")
#load("cem.rda")

summary(all.cem)
if(interactive())
 plot(density(all.cem))

library(cluster)

# diagnostic plots, not used in the paper
if(interactive()){

 plot(density(all.cem),xlim=c(-1000,3500))
 lines(density(all.psm,n=3000),col="red",lty=2)
 lines(density(all.mdm,n=3000),col="blue",lty=3)
 legend(2000,0.0010, legend=c("stratification", "PSM", "MDM"),col=c("black","blue","red"),lty=1:3)

 plot(all.mdm, all.mdm.sd,col="red",xlim=c(-1000,3500),ylim=c(0,3000))
 points(all.psm, all.psm.sd,col="blue")
 points(all.cem, all.cem.sd,col="black")
}


# FIGURE 2 in the paper, file fig3b-col.pdf (Color) in the paper
pdf("fig2-col.pdf",width=11,height=5,pointsize=12)
plot(1/sqrt(all.mdm.n), all.mdm,ylim=c(-1000,4200),xlim=c(0.037,1/sqrt(200)),col=rgb(1,0,0,alpha=0.4),cex=0.7,pch=20,xlab="number matched, scaled as 1/sqrt(matched)",ylab="att estimate",  xaxt = "n")
xx1 <- sort(1/sqrt(unique(c(all.mdm.n,all.psm.n,all.cem.n))), decreasing=TRUE)
xx1 <- xx1[which(xx1<1/sqrt(211))]
axis(1, xx1, (1/xx1)^2)
points(1/sqrt(all.psm.n), all.psm,col=rgb(0,1,0,alpha=0.4),cex=0.7,pch=20)
points(1/sqrt(all.cem.n), all.cem,col=rgb(0,0,1,alpha=0.4),cex=0.7,pch=20)
legend(1/sqrt(225),4300, legend=c("stratification", "PSM", "MDM"),col=c("blue","green","red"),pch=rep(20,3))
e.cem <- ellipsoidhull(cbind(1/sqrt(all.cem.n),all.cem))
polygon(predict(e.cem), col=adjustcolor("gray",0.5))

idx <- which(all.psm.n<800 & all.psm.n>210)
e.psm <- ellipsoidhull(cbind(1/sqrt(all.psm.n),all.psm)[idx,])
lines(predict(e.psm),col="green")

idx <- which(all.mdm.n<800 & all.mdm.n>210)
e.mdm <- ellipsoidhull(cbind(1/sqrt(all.mdm.n),all.mdm)[idx,])
#lines(predict(e.mdm),col="red")
lines(predict(e.mdm),col="red")
dev.off()

# FIGURE 2 in the paper, file fig3b-bw.pdf (Black&White) in the paper
pdf("fig2-bw.pdf",width=11,height=5,pointsize=12)
pch.mdm <- 15
pch.psm <- 8
pch.cem <- 20
col.psm <- rgb(0,0,0,alpha=0.5) #"black" # rgb(0,1,0,alpha=0.5)
col.mdm <- rgb(0,0,0,alpha=0.5) #"black" # rgb(1,0,0,alpha=0.5)
col.cem <- rgb(0,0,0,alpha=0.5) #"black" # rgb(0,0,1,alpha=0.5)

plot(1/sqrt(all.mdm.n), all.mdm,ylim=c(-1000,4200),xlim=c(0.037,1/sqrt(200)),col=col.mdm,cex=0.7,pch=pch.mdm,xlab="number matched, scaled as 1/sqrt(matched)",ylab="att estimate",  xaxt = "n")
xx1 <- sort(1/sqrt(unique(c(all.mdm.n,all.psm.n,all.cem.n))), decreasing=TRUE)
xx1 <- xx1[which(xx1<1/sqrt(211))]
axis(1, xx1, (1/xx1)^2)
points(1/sqrt(all.psm.n), all.psm,col=col.psm,cex=0.7,pch=pch.psm)
points(1/sqrt(all.cem.n), all.cem,col=col.cem,cex=0.7,pch=pch.cem)
legend(1/sqrt(225),4300, legend=c("stratification", "PSM", "MDM"),col="black", #c(col.cem, col.psm, col.mdm),
   pch=c(pch.cem,pch.psm,pch.mdm))
e.cem <- ellipsoidhull(cbind(1/sqrt(all.cem.n),all.cem))
polygon(predict(e.cem), col=adjustcolor("gray",0.3))

idx <- which(all.psm.n<800 & all.psm.n>210)
e.psm <- ellipsoidhull(cbind(1/sqrt(all.psm.n),all.psm)[idx,])
lines(predict(e.psm),col="black",lty=2)

idx <- which(all.mdm.n<800 & all.mdm.n>210)
e.mdm <- ellipsoidhull(cbind(1/sqrt(all.mdm.n),all.mdm)[idx,])
#lines(predict(e.mdm),col="red")
lines(predict(e.mdm),col="black",lty=3)

dev.off()


# FIGURE 2, rescaled a "n" in the x-axis
if(interactive()){
    plot(all.mdm.n, all.mdm,ylim=c(-1800,3500),col="red",cex=.5,xlim=c(-130,722),pch="o",xlab="matched units",ylab="att estimate")
points(all.psm.n, all.psm,col="green",cex=0.5,pch="x")
points(all.cem.n, all.cem,col="blue",cex=0.5,pch="+")
legend(535,3300, legend=c("stratification", "PSM", "MDM"),col=c("blue","green","red"),pch=c("+","x","o"))

e.cem <- ellipsoidhull(cbind(all.cem.n,all.cem))
polygon(predict(e.cem), col=adjustcolor("gray",0.2))

idx <- which(all.psm<3500 & all.psm>-1500)
e.psm <- ellipsoidhull(cbind(all.psm.n,all.psm)[idx,])
lines(predict(e.psm),col="green")

idx <- which(all.mdm<3500 & all.mdm>-1600)
e.mdm <- ellipsoidhull(cbind(all.mdm.n,all.mdm)[idx,])
#lines(predict(e.mdm),col="red")
lines(predict(e.mdm),col="red")
}


# TABLE 1
# The numbers are slightly different from the paper because of randomization of PSM and MDM.
# For stratification no change occurs being the method a deterministic one.

load("psm.rda")
load("cem.rda")
load("mdm.rda")

PSM <- tab.psm[,"est"]
MDM <- tab.mdm[,"est"]
CEM <- tab.cem[,"est"]
tab1 <- rbind( c(min(MDM,na.rm=TRUE), mean(MDM,na.rm=TRUE), quantile(MDM,p=0.5,na.rm=TRUE), max(MDM,na.rm=TRUE)),
c(min(PSM,na.rm=TRUE), mean(PSM,na.rm=TRUE), quantile(PSM,p=0.5,na.rm=TRUE), max(PSM,na.rm=TRUE)),
c(min(CEM,na.rm=TRUE), mean(CEM,na.rm=TRUE), quantile(CEM,p=0.5,na.rm=TRUE), max(CEM,na.rm=TRUE)))
colnames(tab1) <- c("min ATT", "average ATT", "median ATT", "max ATT")
rownames(tab1) <- c("MDM", "PSM", "Stratification")
round(tab1,1)

PSM <- tab.psm[,"sd"]
MDM <- tab.mdm[,"sd"]
CEM <- tab.cem[,"sd"]
tab2 <- rbind( c(min(MDM,na.rm=TRUE), mean(MDM,na.rm=TRUE), quantile(MDM,p=0.5,na.rm=TRUE), max(MDM,na.rm=TRUE)),
c(min(PSM,na.rm=TRUE), mean(PSM,na.rm=TRUE), quantile(PSM,p=0.5,na.rm=TRUE), max(PSM,na.rm=TRUE)),
c(min(CEM,na.rm=TRUE), mean(CEM,na.rm=TRUE), quantile(CEM,p=0.5,na.rm=TRUE), max(CEM,na.rm=TRUE)))
colnames(tab2) <- c("min Std. Error", "average Std. Error", "median Std. Error", "max Std. Error")
rownames(tab2) <- c("MDM", "PSM", "Stratification")

PSM <- tab.psm[,"n"]
MDM <- tab.mdm[,"n"]
CEM <- tab.cem[,"n"]
tab3 <- rbind( c(min(MDM,na.rm=TRUE), mean(MDM,na.rm=TRUE), quantile(MDM,p=0.5,na.rm=TRUE), max(MDM,na.rm=TRUE)),
c(min(PSM,na.rm=TRUE), mean(PSM,na.rm=TRUE), quantile(PSM,p=0.5,na.rm=TRUE), max(PSM,na.rm=TRUE)),
c(min(CEM,na.rm=TRUE), mean(CEM,na.rm=TRUE), quantile(CEM,p=0.5,na.rm=TRUE), max(CEM,na.rm=TRUE)))
colnames(tab3) <- c("min n", "average n", "median n", "max n")
rownames(tab3) <- c("MDM", "PSM", "Stratification")

library(xtable)

xtab1 <- xtable(round(tab1,1),digits=1)
xtab2 <- xtable(round(tab2,1),digits=1)
xtab3 <- xtable(round(tab3,0),digits=0)

t1 <- capture.output(print(xtab1, floating=FALSE, include.rownames=FALSE))
t2 <- capture.output(print(xtab2, floating=FALSE, include.rownames=FALSE))
t3 <- capture.output(print(xtab3, floating=FALSE, include.rownames=FALSE))
t1 <- c(t1[-(1:2)],"\n")
t2 <- c(t2[-(1:2)],"\n")
t3 <- c(t3[-(1:2)],"\n")
ttt <- c("\\documentclass[11pt, oneside]{article}\n","\\begin{document}\n",
"\\begin{table}\n", "\\begin{tabular}{c}\n",
t1, "\\\\\\\\",t2,"\\\\\\\\", t3,"\\end{tabular}\\caption{Distribution of estimated ATTs, their standard error and number of matched units for the different matching methods.}\\end{table}\n\\end{document}")
myfile <- file("table1.tex","w")
cat(ttt,file=myfile)


